Spin liquid ground state in a two dimensional non-frustrated spin model 
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We consider an exchange model describing two isotropic spin-1/2 Heisenberg antiferromagnets 
coupled by a quartic term on the square lattice. The model is relevant for systems with orbital 
degeneracy and strong electron- vibron coupling in the large Hubbard repulsion limit, and is known 
to show a spin-Peierls-like dimerization in one dimension. In two dimensions we calculate energy 
gaps, susceptibilities, and correlation functions with a Green's Function Monte Carlo. We find a 
finite spin gap and no evidence of any kind of order. We conclude that the ground state is, most 
likely, a spin liquid of resonating valence bonds. 

PACS Numbers: 75.10.Jm, 75.40.Mg, 75.50.Ee, 02.70.Lq 
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The existence of a homogeneous resonating valence 
bond (RVB) ground state (GS) for two-dimensional (2D) 
spin-1/2 systems has been a subject of intense theoretical 
study in the past two decades. The idea that the reduced 
dimensionality and the small spin value might enhance 
quantum fluctuations up to the point of destroying the 
classical Neel antiferromagnetic order was first put for- 
ward by Anderson in 1973. IT] 

The subject became a hot topic after Anderson sug- 
gested that the parent compounds of the high-Tc copper- 
oxides - argued to be well described by a S— 1/2 Heisen- 
berg model on a square lattice - might have a spin liquid 
GS, and that superconductivity would result from doping 
such a spin liquid. 

Since then, many studies on the S=l/2 2D Heisenberg 
model have shown that reduced dimensionality and low 
spin are not sufficient to stabilize a RVB GS. First, it is 
now well established that the square lattice case is Neel 
ordered, [pj More surprisingly, three-sublattice Neel order 
is also likely to survive on the triangular lattice, where 
the model is frustrated. 0,0] A stronger lattice frustration 
- like in the kagome lattice g) -, multiple-spin exchange 
terms on the triangular lattice Ig], or frustration due to 
the coupling constants - like in the Ji — J2 Heisenberg 
model on the square lattice |7|j^ - might be more effective 
in stabilizing a spin liquid GS. The evidences for a triplet 
gap and the definite assessment about the liquid nature 
of the GS are, however, either limited to small lattice 
exact diagonalizations (iV = 36 sites) , ^^ or must cope, 
when Monte Carlo is used, with the sign problem. ||] 

We propose here a clean example of a RVB ground 
state realized, in absence of any frustration, in a 2D 
model where an extra spin-1/2 T, representing an orbital 
degree of freedom, is coupled to the usual spin S. Specif- 
ically, the model we have considered is a spin-exchange 
Hamiltonian of the following form: 



where J > 0, the summation runs over the nearest- 



neighbor (nn) sites of a square lattice, X. 
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1/2), and Xiy - (2T, 1 
forth referred to as ST model, describes for J > the 
low-energy physics of an insulating crystal with one elec- 
tron per site in a two-fold degenerate orbital, in the 
limit of large on-site repulsion (Mott insulator) and in 
presence of Jahn- Teller (JT) effect. [^ The derivation 
of Eq. (|l|), in the same spirit of the derivation of the 
Heisenberg model from the large-U Hubbard limit, is 
standard. [|l^,^ 11 1 The crucial physical condition to be 
verified is that, among the possible two-particle states 
obtained upon virtual hopping, the inter-orbital singlet 
should turn out to be the lowest in energy, which is in- 
deed the case when a strong dynamical JT effect is at 
play. A different, perhaps more common, physical situa- 
tion is that the lowest two-particle intermediate state is 
a Hund's triplet, in which case the exchange model leads 
to spin ferromagnetism. P^jl^l In the general case, the 
exchange Hamiltonian contains pseudo-spin anisotropic 
terms, and some aspects of its phase diagram have been 
addressed with different techniques in Refs. [|lO|Jl3| ]. 

A class of SU(n)-invariant generalizations of the 
Heisenberg model has been previously proposed as candi- 
date non-Neel antiferromagnets. ||lj,|l5[ It turns out that 
the ST model is unitarily equivalent to a particular SU(4) 
model in that class, with Uc = m = 1 (notations of Ref. 
||l5|). The n ^ cxD limit of these SU(n) models has been 
studied in detail, and shown to have, for Uc — m = 1, 
a spin-Peierls GS both in 2D, ||l5|-|l7|| as well as in one 
dimension (ID). [Q Recently, we have verified that the 
spin-Peierls nature of the GS persists, in ID, down to 
the ST model point, i.e., n = 4. |H| The behavior in 2D is 
more subtle. We show in the present letter that a RVB 
GS is found for n = 4 in 2D, in marked contrast with the 
n — > 00 limit predictions. 

The relationship between the ST model and the usual 
Heisenberg model - i^Hcis = (Jf^) J^Ui) -^i j ~ ^^ made 



more clear by working in the singlet subspace with the 
overcomplete basis set of the valence bond (VB) config- 
urations, i.e., products of non-overlapping singlet dimers 
(i, j) connecting pairs of sites on opposite sublattices and 
covering the lattice. More precisely, if {i,j)s(T) = (TUj 
— iiTj)s(T) denotes a singlet bond between sites i and 
j for the S (T) variables, we restrict the VB configura- 
tions considered to the invariant subspace - containing 
the GS, by positivity arguments - in which every sin- 
glet dimer {i,j) is a product of an S- and a T-singlet, 
Ihj) = {hJ)s{i,j)T- Since (i,j)s(T) = -(j,J)s(T), we 
assume the sign convention that the leftmost index in the 
singlet pair always belongs to sublattice A. It is possible 
to show that the GS of both the ST and the Heisenberg 
model can be written as linear combinations with posi- 
tive coefficients of VB configurations, with the previous 
restrictions and sign conventions. The construction of ex- 
cellent variational wavefunctions based on combinations 
of VB configurations has been shown to be possible even 
for the Heisenberg case provided sufficiently long-ranged 
bonds are allowed. |18| On the other hand, it is quite 
clear that whenever short-ranged bonds are the domi- 
nant ones, the GS will have a characteristic length ^ and 
there will be a finite gap, in the thermodynamic limit, to 
the lowest triplet excitations. 

The bond operators x\f = (28^ • Sj - 1/2), and 

(T) 

similarly X^ ■ , have simple properties when acting on 
VB configurations, because they affect at most two sin- 
glet pairs. Indeed, it is known that Q] Xl-'{i,j)s — 

-{k,j)s{i,l)s- Us- 



-2{i,j)s, andxf^)(i,j)5(fc,/)s 



ing these rules for X^ , and X^ , , it is straightforward 

(S) (T) 

to show that the bond operators Xi,j ~ —X^.X^, ap- 



pearing in the ST model, Eq. (hi) 
rules: 



obey the following 



X3,kiiJ){k,l) = -{k,j){i,l) , 



(2) 



with n = 4. Notice that these relationships are formally 
identical to those relevant to the Heisenberg case, except 
for a coefficient n = 4, in place of n = 2, when a nn 
Hamiltonian bond (ij) acts on a single dimer (i, j). This 
enhanced coefficient favors the formation of short-ranged 
bonds in the GS of the ST model, making the suppression 
of Neel long range order (LRO) more likely. In the limit 
n — > oo we recover, from Eq. (0), the known results both 
in ID, WM and 2D. fia] Indeed, for n ^ oo the only sur- 
viving VB configurations are those with nn dimers only. 
In ID this leads to a doubly-degenerate spin-Peicrls GS. 
WM In 2D the model maps |l5[ onto the purely kinetic 
limit of the quantum dimer model (QDM) of Rokhsar and 
Kivelson. |19J (This purely kinetic point of the QDM is 
believed to be characterized by a plaquette resonating 
VB state, [O] breaking translational invariance.) 

By working with the VB basis, the action of the 



Hamiltonian Hipi on any basis element ipi defines a non- 
symmetric matrix, Hipi = ^ hj^iipj, with all elements 
hj,i non-positive, as implied by Eq. (|^). The right eigen- 
vector of hj,i with minimum eigenvalue, corresponding 
for hj^i < to the GS of i?, can be computed by ap- 
plying the power method, as implemented stochastically 
by means of the Green's Function Monte Carlo (GFMC) 
method. The GFMC is in fact not limited to symmetric 
matrices, and there is no sign problem when all hj^i are 
non-positive. B Using the VB basis, the GFMC turns 
out to have extremely small statistical errors, compared 
to the more conventional algorithm [Q employing an 
Ising basis. In this formulation, the GFMC does not re- 
quire the calculation of the overlaps (V'i | 4'i) between 
VB configurations. Details of the method are given else- 
where. |2^] 

This new and simple GFMC allows us to obtain a very 
accurate determination of the triplet gap by performing 
two independent simulations of the singlet GS and the 
triplet lowest excited state. In the latter case, the basis 
employed is slightly modified, by allowing for the pres- 
ence of a single triplet bond (i,j)* = (tiij + ii]j)s{h3)T 
in each element of the VB basis. The corresponding 
rules for the application oi Xi^j are: Xi,j{i,jy — 0, and 
Xj^k{hjY{k,l) = —{k,j){i,lY. Notice that this implies 
the absence of sign problem in the triplet subspace as 
well. 

Fig. n shows the results obtained for the triplet gap 
for L X L square lattices with L up to 24, as a func- 
tion of the inverse volume. The corresponding data 
for the Heisenberg case, [EOJ also shown for compari- 
son, are consistent with a finite size gap A^ scaling to 
zero as a/L'^ + b/L^ + ■ • •. The dashed lines (see also 
left inset in Fig. 1^) show our best two-parameter fit to 
the ST data obtained by imposing the same gap be- 
havior as in the Heisenberg case: such a fit is clearly 
unsatisfactory. Instead, the solid line through the ST 
data is the result of a three parameter fit of the form 
Al = a + a/L^ + b/L^ + • • •, [^ giving a clear evidence 
of a finite gap A in the thermodynamic limit. Clearly, 
the detailed finite-size behavior of the ST gap in 2D is 
non trivial, requiring the simulation of quite large lat- 
tices {L = 24) to pin-down the presence of a gap. This 
suggests the presence of a length scale ^ of the order of 
10 -^ 20 lattice spacings. This behavior should be con- 
trasted to that of the triplet gap for the ST model in ID, 
shown in the right inset of Fig. ^, where the size scaling 
of the gap is straightforward {£, « 1). 

In principle, either a VB crystal with some broken 
spatial symmetry, as in ID, or a homogeneous spin liq- 
uid is compatible with the existence of a spin gap. In 
order to investigate the possible kinds of LRO which 
might characterize the GS of the ST model in 2D, we 
have calculated the expectation value of several spin-spin 
correlation functions by the forward walking technique. 
p3|,p0[ Fig. shows the results for the spin structure fac- 



tor 5(77, tt) = X]r(~l)''(^o • Sr), which should diverge 
with the system vohime L^, when the Neel order pa- 
rameter ra^ is non- vanishing, S'(7r, tt) ^ [Lm'^Y' . The 
comparison between the Heisenberg and the ST model 
results for S{'n, tt) strongly suggests the absence of Neel 
order in the ST case, in agreement with the presence 
of a triplet gap. More interestingly, the inset in Fig. 
shows the results obtained, in the ST case, for the dimer- 
dimer structure factors with nn bonds in the aj-direction, 
S'^-^d) = J:^ e^^-^S^S^+,S^S^+,). The q = (^, tt) and 
q = (0, tt) dimer structure factors are clearly finite in the 
L ^ 00 limit; the q — (tt, 0) also does not seem to di- 
verge linearly with the volume, but the present data are 
limited to too small sizes {L < 10) to be conclusive. 

In order to have further evidence about the existence 
of a true homogeneous spin liquid GS, we have directly 
calculated, using the new more efficient GFMC algo- 
rithm described previously, the response of the system 
to symmetry-breaking operators. It is practically im- 
possible to exclude all possible types of crystalline or- 
der numerically; in the following we restrict our consid- 
eration to the most plausible types of order, involving 
either a broken translation symmetry T with momenta 
compatible with a real GS (q = (7r,0), (0,7r), (7r,7r)), 
or/and a broken tt/2 rotation symmetry R. This in- 
cludes all types of crystalline dimer and plaquette or- 
der proposed so far. p9 16 1^ More precisely, we have 



perturbed the ST Hamiltonian by adding a term aO, 
with O an operator which breaks one of the symmetries 
above: O — J2r ^^'^'^^r,r+x with the appropriate q for 



the translation case, and O — J^ri-^J^.r 



X, 



r,r+y 



) for 



the tt/2 rotation. On finite size, the GS expectation value 
of O vanishes by symmetry, and the GS energy per site 
has corrections proportional to a^, Eq = eo — Xoa^/2, 
Xo being the generalized susceptibility associated to the 
symmetry-breaking operator O. On the other hand, if 
symmetry breaking occurs in the thermodynamic limit, 
it is possible to have limc^o limi^oo(0/L^) = p ^ 0. 
In the latter case, xo has to diverge as i ^ cxd. More 
precisely, it is possible to show that xo is bounded from 
below by the order parameter times the system volume 
squared, xo > const p^L^. These susceptibilities are 
therefore a very sensitive tool - much more than the 
usual square of the order parameter ~- for detecting LRO. 
For instance, as shown in Fig. for the ST model in 
ID, the presence of dimerization in the thermodynamic 
limit ||9(] is readily inferred from the behavior of xo , with 
O — X]i(~l)*^i,i+ii even for very small system sizes 
{L < 12). Fig. ^ shows the results obtained for the ST 
model in 2D. For all the symmetry breaking operators 
considered, the associated susceptibilities are found to be 
finite. The largest susceptibility is found to be the one 
associated to columnar order (T(7r, 0) in Fig. Q), and is 
only weakly increasing with size, eventually saturating to 
a constant, in marked contrast to the strong divergence 



of the ID analog (see Fig. ||). 

In conclusion, we have studied a non-frustrated ex- 
change Hamiltonian, Eq. (|l|), which describes the low- 
energy physics of a Mott insulator with orbital degen- 
eracy in the regime in which the inter-orbital singlet is 
the lowest-energy intermediate state available to virtual 
hopping. We find a clear evidence for a spin gap with- 
out crystalline VB order in 2D. A homogeneous liquid of 
resonating valence bonds appears as the natural candi- 
date GS for this model. An important question about 
the nature of the low-lying excitations remains to be ad- 
dressed: is there a branch of gapless excitations in the 
singlet sector, or the system is gapped to all excitations? 
The answer to this question will require future work. 
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FIG. 1. Finite size gap between the GS and the first 
excited triplet state for the ST model and for the Heisenberg 
model |2(| on the square lattice. The ST data are obtained 
by the new GFMC method working with the VB basis. The 
solid and dashed lines are different fits described in the text. 
The data for the ST model clearly indicate the existence of 
a spin gap in the thermodynamic limit. Left inset: The ST 
data plotted versus 1/L. Right inset: The triplet gap for the 
ST in ID. 
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FIG. 3. The susceptibility to breaking translation sym- 
metry for the ST model in ID, which is known to have a 
spin-Peierls-Iike dimerization. Data are obtained by exact di- 
agonalizations. Inset: the extrapolated a = value of xo 
divided by the square of the system volume L^ . 
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FIG. 2. Finite size structure factors obtained by GFMC 
with the forward walking technique. From the spin-spin 
S{tv,tt) we find no clear sign of Neel LRO. The size scaling 
of the Neel order parameter for the Heisenberg model [M is 
shown for comparison. Inset: The important components of 
the Fourier transform of the dimer-dimer correlation functions 
appear to be finite. 



FIG. 4. Susceptibilities to breaking the most important 
crystal symmetries for the ST in 2D. The finite-size, finite-a 
data for xo can be readily extrapolated, by quadratic fits, to 
a = values: notice the distinctly different behavior for the 
corresponding ID results of Fig. 0. The insets summarize the 
results obtained for the extrapolated xo- 



